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Abstract 

We investigated the dynamical evolution of clusters of galaxies in virial equilibrium using Fokker-Planck 
models and self-consistent A^-body models. In particular, we focused on the growth of a common halo, 
which is a cluster-wide halo formed by matter stripped from galaxies, and the development of a central 
density cusp. The Fokker-Planck models include the effects of two-body gravitational encounters both 
between galaxies and between galaxies and common halo particles. The effects of tidal mass stripping from 
the galaxies due to close galaxy-galaxy encounters and accompanying dissipation of the orbital kinetic 
energies of the galaxies were also taken into account in the Fokker-Planck models. We find that the results 
of the Fokker-Planck models are in excellent agreement with those of the A^-body models regarding the 
growth of the common halo mass and the evolution of the cluster density profiles. In the central region 
of the cluster, a shallow density cusp, approximated by p{r) cx r~" {a ~ 1), develops. This shallow cusp 
results from the combined effects of two-body relaxation and tidal stripping. The cusp steepness, a, weakly 
depends on the relative importance of the tidal stripping. When the effect of stripping is important, the 
central velocity dispersion decreases as the central density increases and, consequently, a shallow {a < 2) 
cusp is formed. In the limit of no stripping, usual gravothermal core collapse occurs, i.e. the central 
velocity dispersion increases as the central density increases with a steep {a > 2) cusp left. We conclude 
from our consideration of the origin of the cusp demonstrated here that shallow cusps should develop in 
real galaxy clusters. 
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dynamics 



1. Introduction 

Many observations suggest that clusters of galaxies contain substantial amounts of mass in the form of a background 
intergalactic stellar component (e.g., Vflchez-Gomez et al. 1994; Scheick, Kuhn 1994; Ferguson et al. 1998 and references 
therein). It has usually been thought that such a component has its origin in stars stripped from galaxies by galaxy- 
galaxy encounters, or by the cluster tidal field. However, detailed numerical simulations of galaxy clusters are required 
to confirm how efficient such mass stripping is and how the structure of the background component evolves. 

It was only recently that fully self-consistent A^-body simulation of galaxy clusters, with each galaxy represented 
by many particles, was made possible, because such a full simulation requires huge computing resources. Funato et 
al. (1993) investigated the evolution of isolated clusters of galaxies using self-consistent A^-body simulations. Similar 
simulations were also performed by Bode et al. (1994) and Garijo et al. (1997). More recently, Sensui et al. (1999, 
SFM99; 2000, SFMOO) followed up the study of Funato et al. (1993) by performing larger A^ simulations for a wider 
variety of initial conditions. In their simulations, clusters were initially in virial equilibrium, and all of the cluster mass 
was initially bound to individual galaxies; a cluster was composed of about 100 identical galaxies, and each galaxy was 
modeled with 500-4000 particles. SFM99 and SFMOO found that, in all cases, more than half of of the total mass is 
stripped from galaxies in a few crossing times of the cluster, and that the stripped matter forms a smooth cluster-wide 
common halo. They also found that a density cusp, approximated as pcx r~°' {a being form —1 to —1.5), develops in 
the central region. However, the physical mechanism of the cusp formation was not clear. 

Even if we are allowed to regard clusters of galaxies as "pure" self-gravitating many-body systems composed of 
point masses (stars and dark matter particles), we have to admit that their evolution is complicated by many phys- 
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ical processes which take place simultaneously: gravitational encounters, physical collisions, and merging between 
galaxies; galaxy-background particle interactions, mass stripping from galaxies, internal evolution of galaxies, etc. 
Therefore, there is no doubt that fully self-consistent A^-body simulations arc most desirable to obtain qualitatively 
and quantitatively reliable results on cluster evolution (see Garcia-Gomez et al. 1996). 

Nevertheless, more idealized or more simplified models that include (part of) relevant physical processes in approx- 
imate ways, such as Fokker-Planck (FP, hereafter) models which we used in the present study, are still very useful 
regarding two points. First, we can quickly obtain an overview of cluster evolution based on various initial conditions 
by using approximate models, which are inexpensive in computation time. Second, it is expected that we can un- 
derstand more clearly the role and importance of each physical process, because we can easily separate the effects of 
various processes in approximate models. On the other hand, it is often difficult to "understand" the results of iV-body 
simulations, where all processes are automatically included and connected through a complicated procedure. Thus, 
iV-body and approximate models arc complementary to each other, and comparison between these models will be 
particularly helpful to our further understanding of the physics of cluster evolution, regardless of whether agreement 
or disagreement is found. 

Merritt (1983, 1984, 1985) explored the evolution of clusters of galaxies using such approximate models based on 
an isotropic orbit-averaged FP equation (Cohn 1980). The FP equation describes the evolution of the phase-space 
distribution function driven by gravitational two-body encounters, and has been applied mainly to globular clusters 
(see, e.g., Spitzer 1987 for a review). A typical cluster of galaxies is composed of iV = 100-1000 galaxies, and thus 
the time scale of two-body relaxation between galaxies, Tr, is of the same order as the crossing time of the cluster, 
Tcr (Tr/Tcr ~ O.liV/luiV). Thcrcforc, two-body relaxation is expected to play an important role in the evolution of 
galaxy clusters. In his models, Merritt also took account of mass stripping from galaxies due to tidal interactions 
between galaxies. Since the stripped particles (stars or dark matter particles) were considered to be much lighter than 
galaxies, these particles were treated as being massless. Dynamical friction between the galaxies and the background 
massless particles was also included. Thus, in Merritt's models the evolution of the galaxies and the background was 
treated self-consistcntly. FP models arc continuum models, and therefore have an advantage that they are free from 
statistical fluctuations, which arc unavoidable in A^-body models and make it difficult to understand the underlying 
physics. Yepes et al. (1991) and Yepes and Dommguez-Tenreiro (1992) studied the evolution of galaxy clusters by 
solving the moment equations of the FP equation. They also included a background component of massless particles 
in their models, but did not take account of tidal stripping. 

The purpose of the present work was to investigate the dynamical evolution of clusters of galaxies after virialization. 
In particular, we aimed to clarify the physical mechanisms of the growth of the common halo, and of the development 
of the central density cusp, which were investigated by SFM99 and SFMOO with A^-body models. For this purpose, we 
used FP models, which arc based on Merritt's (1983) formulation, but arc improved in many aspects, and compare 
them with the self-consistent A^-body models of SFM99 and SFMOO. As shown below, the results of the FP and A^- 
body models agree with each other surprisingly well. This suggests that our FP models correctly include all important 
physical processes, and thus can be used as a useful tool for studying the dynamics of galaxy clusters. 

This paper is organized as follows. The basic equations of the FP models are given in section 2. Section 3 describes 
the initial conditions of our simulations, and section 4 briefly describes the numerical schemes used for our FP and 
A'-body simulations. The simulation results are given in sections 5 and 6; in section 5, the FP and A^-body models 
are compared with respect to the growth of the common halo and the development of the central cusp; in section 6, 
the physical mechanism of the cusp formation is examined in detail based on FP models. We summarize and discuss 
our results in section 7. 

2. Fokker Planck Models 

In this section, we describe FP models of clusters of galaxies. In FP models, individual galaxies are treated as 
point particles, except that the effects of tidal stripping are considered. Detailed FP simulations of galaxy clusters 
were first performed by Merritt (1983). We closely follow his formulation. There are, however, three major differences 
between his and our models: (1) Merritt used isotropic FP models, but we used anisotropic FP models. (2) We 
used a new cross-section formula of mass stripping from galaxies based on recent TV-body experiments by Funato and 
Makino (1999, FM). (3) Dissipation of the orbital kinetic energies of galaxies during galaxy-galaxy encounters was 
also considered, while it was neglected in Merritt's simulations. 

2.1. The Fokker-Planck Equation 

Gonsider a self-gravitating many-body system composed of different mass components with particle mass rui (i = 
0,1,2,...). We write the number density in ^-space, i.e. the distribution function, of component as fi(r,v,t). By 
assuming that the system is spherically symmetric and in dynamical equilibrium, fi can be written as a function of 
the energy per unit mass, E, and the angular momentum per unit mass, J. Then, the orbit-averaged FP equation, 
under the fixed gravitational potential, can be written as follows (Gohn 1979; Takahashi 1997): 
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+ ^(Dmf,+DnE^% + D^^,% ] , (1) 



where R denotes the scaled angular momentum, A is a weight function, and Dei, Deei, etc. are flux coefficients. See 
Takahashi (1997) for a complete description of this equation. To follow the evolution of the system, this FP equation 
is solved in combination with Poisson's equation (Cohn 1979). 
We introduce the mass density in phase space, 

9i = inJ.,. (2) 

Here, we choose the i = component for representing the background, or a common halo. We can assume that the mass 
of a background particle (star) is negligible compared to that of a galaxy. However, the mass density of background 
particles is comparable to that of galaxies. We now take the limit mp — > while keeping go finite (i.e. /o ^ oo). In 
the FP equation (1), component mass rrii appears only in the flux coefficients. In the limit mo 0, the expressions 
for the first and second-order flux coefficients look like 

DE(i = 0, DEiCcmi^gj (i = 1,2,3, .. .), (3) 
Deez, DERi, DREi, DRm(xy^mjgj (z = 0, 1, 2, . . .). (4) 

Note that we write here only the dependence on to, and gi symbolically. Exact expressions can be easily obtained 
from the equations given in the appendix of Takahashi (1997). 
The FP equation for the background component {i = 0) becomes 

.dga d ( ^go ^5oV d , ^ dgo\ 

^^^dE [^''''"dE + M ) + M [^'"'^'dE + M ) ' 

which is obtained by multiplying equation (1) by toq E^nd then taking the limit niQ — > 0. Since the background particles 
are massless, they are completely coUisionless, i.e., they do not interact with each other. However, galaxies suffer 
dynamical friction from the background particles that have a finite mass density. Thus, the orbital energy of the 
galaxies is transferred to the background. 

2.2. Mass Stripping Rates 

Consider an encounter between a test galaxy of mass m and a field galaxy of mass m/, with impact parameter p 
and relative velocity V. FM showed, both numerically and analytically, that the relative change of the test galaxy's 
mass in an encounter can be approximated by 

TO \ P / ^ V / \ m J 

with rjr ~ 2-3 and ry^ ~ 2-3, for the cases with p^r^, y ^ Ug, and m = mf. Here, C is a dimensionless constant, and 
Tg and iig are the size (represented by the virial radius in the following) and the three-dimensional internal velocity 
dispersion of the test galaxy, respectively. Since FM did not carry out simulations of encounters between unequal-mass 
galaxies, the value of rjm is not determined from their simulations. In the impulse approximation, the velocity change 
of a star. Aw, is proportional to ruf. Therefore, if the effect of the term (Aw)^ is dominant in removing mass from the 
test galaxy, we may expect rjm — 2; on the other hand, if the effect of the term vAv is dominant, rfm — 1 (sec Mcrritt 
1983; FM). 

Now suppose that a test galaxy experiences successive encounters with field galaxies of mass ruf. The rate of change 
in the test galaxy's mass is given by 

where nf is the number density of the field galaxies and F{V) is the (normalized) distribution function of the relative 
velocity between the test and field galaxies; Pmin and Pmax are the minimum and maximum impact parameters, 
respectively. If we assume a Maxwellian distribution with a three-dimensional velocity dispersion, V/, for F, this 
equation becomes 



^dt \3J ^'"/"g^/ vf ) / x'-^^dxj y'-^'^eM-yV2)dy. (8) 
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Here, we define x = p/r^ and y = ^/3V/Vf. The first integral becomes ln(a::niax/2^min) for ijr = 2, and l/xmin — l/a;max 
for rjr = 3. Therefore, if we take Pmin and Pmax as the size of the galaxy and that of the cluster, we may consider that 
the integral remains roughly constant during the cluster evolution for rjr = 2-3. The second integral is also a constant 
of the order of unity for rj^ = 2-3; it equals 1 when — 2, and (7r/2)^/^ when rjy = 3. 
Therefore, we may rewrite equation (8) as 

where C is a dimensionless constant factor and pf = nifUf is the mass density of the field galaxies. Note that the 
definition of C depends on rjr and ijv When rjv ^ 2, for example, C = (87r/3)^/^ ln(pniax/pmin)C" for rjr ~ 2, and 
C = (87r/3)i/2(l/p„i„ - l/pniax)C" for rjr = 3. 

Equation (9) includes two parameters of galaxy's internal structure, Tg and v^. Because FP models do not follow 
the internal evolution of galaxies, we have to adopt some relations to express the changes in these parameters by m 
only. One such relation should be the virial theorem, w| oc Gm/r^. As the other relation we assume Wg (x with C ~ 
1/4-1/3 (see FM and SFM99). Using these relations we may write the two parameters as 



Wg / TO \ ^ Tg / TO ^ ^ 

Wgo Wgo/ ' ''go V'7^gO, 



(10) 



Finally, we obtain 
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where a sum over all galaxy components / is taken. We then orbit-average (see, e.g., Takahashi 1997) equation (11) 
to obtain the mass-loss rate as a function of E and J. 

2.3. Energy Dissipation Rates 

Stars escape from galaxies via galaxy-galaxy encounters. This is possible because part of the orbital kinetic energies 
of galaxies is converted into the energy of internal motions of the stars in the galaxies. That is, these encounters are 
inelastic. If the sum of the galaxy binding energies is much smaller than the cluster binding energy, i.e., if (wg/V^ci)^ <C 1 
(V^i is the cluster velocity dispersion), the effect of energy dissipation due to inelastic encounters is not very important. 
However, for our model clusters shown below, this condition is not well satisfied [(iig/Vci)^ = 5/32], and therefore energy 
dissipation is expected to play a non-negligible role. 

The binding energy of a galaxy (in virial equilibrium) is £g = mv'i/2. The rate of its change is therefore given by 



dhiea- d\nm d\nv ^dhi 



dt dt ' dt ^ ^ dt ^ 

Here, in the second equality, we used equation (10). The mass-loss rate is given by equation (11). 

The decrease in the binding energy of a galaxy should be compensated by a decrease in the orbital kinetic energies 
of the galaxies involved in encounters with that galaxy. For simplicity, we assume that the change in the binding 
energy of a galaxy is just equal to the change in its orbital energy, and that the stripped stars have the same orbital 
energy per unit mass as that of the parent galaxy. This is valid at least on average. The dissipation rate of the orbital 
energy of galaxy to is given by 

dt J dt m dt 

(remember E represents the energy per units mass). The above assumption ensures that the total energy of a galaxy 
cluster, including the internal energies of the galaxies, is conserved. 

2.4. Implementation of Mass Stripping and Energy Dissipation 

The distribution function evolves owing to mass stripping from galaxies as well as owing to two-body relaxation. 
The effect of mass stripping can be included in the FP equation by adding a source (or sink) term (see, e.g., Quinlan, 
Shapiro 1989). Thus, in this case the FP equation becomes 

A^ = T, + S,, (14) 

where F,; denotes the FP collision term [the right-hand side of equation (1)] and Si represents the source term due 
to mass stripping. Wc neglect merging events between galaxies for simplicity. For the initial conditions which we 
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adopted in this work, TV-body simulations showed that mergers were rare (SFM99; SFMOO). Therefore, neglect of 
merging events would not cause serious errors for our models. The source term is calculated in practice as follows. 

For each galaxy component i and for each grid point in {E,R) space, we can calculate the change in mass Ami{E , R) 
during a time step At from the orbit-averaged stripping rate, {dm / dt) qa- Thus, we know that galaxies of mass m ~ mi 
at time t and position {E, R) will have mass m' = + Am^ (< m) at time t + At, and that each of these galaxies 
provides mass m — m' for the background. The number of the galaxies does not change during this mass loss, unless 
they are completely disrupted. 

Since we use discrete mass components, a new mass m' is generally not equal to any components mj, and therefore 
we have to somehow distribute the galaxies of mass m' among some components. As Lee (1987) and Quinlan and 
Shapiro (1989) did for stellar mergers in star clusters, we distribute the galaxies by linear interpolation: if n galaxies 
of mass m' are to be formed and if mj < m' < Wj+i, then n{mj^i — m')/(mj+i — mj) galaxies are added to the 
mj component and n{m' — mj)/ (wj+i — mj) galaxies are added to the m^+i component. If m' < mi, where mi is a 
minimum galaxy mass assumed, the galaxies are distributed between mp = (background) and mi. Further we assume 
that the galaxies after mass stripping and the ejected background particles have the same orbital energy and angular 
momentum per unit mass {E, R) as the galaxies had before stripping. The effect of energy dissipation is considered 
separately, as described below. 

By doing the above procedure over all of the galaxy components at each grid point in {E,R) space, we can calculate 
the change in the distribution function due to mass stripping, i.e., the source term. Si. Note that this procedure 
guarantees conservation of the total mass and energy of the cluster. 

Merritt (1983) formulated mass stripping as a first-order "mass diffusion" term in the FP equation [see his equation 
(22)]. In his formulation he assumed a distribution function, f{E,m,t), which is continuous with respect to m as well 
as E. In integrating the FP equation numerically, he used, of course, discrete mass grids. On the other hand, we have 
assumed discrete mass components from the beginning. Although Merritt's formulation and our formulation may look 
different apparently, they both consider only the first-order term (Am) and are essentially the same. 

The effect of energy dissipation (cooling) is included in the FP equation in the same way as the effect of binary 
heating in star clusters (see, e.g., Takahashi 1997). That is, the cooling term, which may be calculated from equation 
(13), is added to the first-order energy coefficient, De- In practice, the amount of energy dissipation during At is 
calculated so that it should be consistent with the above-described way of distributing stripped galaxies among mass 
components. 

2.5. Time Scales 

A commonly used definition of the two-body relaxation time is given by Spitzer (1987, equation [2-61]) as tr = 
(7j^y3)/((Aw||)^)i,=i,j^. Here, a Maxwellian velocity distribution of field particles is assumed and is the three- 
dimensional velocity dispersion; {(Av\\)'^)v^u^ is the mean square value of the change per unit time in the velocity 
component parallel to the initial velocity of a test particle, evaluated at v ~ Vm- For the present cases, assuming that 
all galaxy components follow Maxwellian velocity distributions with the same velocity dispersion, Vd, we obtain 

_ _ 0.065 Vj ^^^^ 

G'HnA^mfpf inA (m/)pPgai ' 

/ 

where pgai ~ "^Pf is the total density of the galaxies and {mf)p ~ '^'^fPf/Pgai is the mass-weighted mean galaxy 
mass. Equation (15) represents the time scale of the velocity diffusion due to encounters between galaxies. Similarly, 
we define the dynamical friction time as idf = |'Um/(Awj[)t,=„,^ |, where (Av||)„=„^^ is the mean change per unit time 
in the velocity component parallel to the initial velocity, evaluated at v ~ Wm. On the same assumptions as used in 
deriving equation (15), we find that the dynamical friction time for a galaxy of mass m is given by 

O.UVj 0.13 Vj 

G2 1nA^(m + m/)p/ G^lnA (m/)pPgai + mptot 
/ 

where ptot is the total density including both the galaxy and halo components. Unlike the relaxation time (15), the 
dynamical friction time (16) depends on the test-particle mass and on the halo density. From equations (15) and (16), 
it follows that 

0.5(1 + ^^). (17) 

In single-component clusters, = T^f. In multicomponent clusters, Tdf can be much smaller than for massive 
galaxies; this holds also for light galaxies, if ptot 3> Pgai- 
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We define the tidal mass stripping time as tts = \m/{dm/dt)\. Using equation (11) with set equal to 2V^^, we 
find 



- ^ 2 „,V. / ,,„,-l\ ■ (18) 

(TO/ )pPgal 

Using equations (16) and (18), with rgo = G'mgo/2UgQ, we obtain 

71^ ^ 311nA J77_y''~"^'"^ / WgoV""" (TO/)pPgal+TOptot . ^ 

For our standard set of parameters, rjm = 2, ?/„ 2, and = 1/4 (see section 5), it follows that 



This indicates that Tts ~ Tdf in our initial models (see section 3). As the fraction of the common halo increases, tidal 
stripping becomes less important (Merritt 1983). 



3. Initial Conditions 



We use the same initial cluster model as that in SFM99. That is, the initial model is a Plummer model composed of 
A"g = 128 identical galaxies with no common halo. Exceptions are models FPAl and FPA2, which we describe later. 
The initial galaxy model is also given by a Plummer model as in SFM99. The ratio of the virial radius of each galaxy 
Tvr to the virial radius of the cluster i?vr is set equal to 1/20. 

Note that initial galaxy models affect FP models only through the mass-stripping rate, equation (11). This equation 
includes two parameters of the internal structure of galaxies, rgo and Wgo, but these are related by the virial theorem 
(we take Tg = Tvr) for a given galaxy's mass. Hence, in FP models the structure of galaxies is represented by only one 
parameter, e.g., Tvr. In reality, even if the virial radii are the same, different galaxy models will give different stripping 
rates. Such differences in galaxy models could be reflected in FP models only by adjusting an undetermined numerical 
constant, C, in equation (11). 

We use a system of units in which G = 1, mg = 1, and £g = 1/4, where G is the gravitational constant (Heggie, 
Mathieu 1986); hence r^^ = Gm^/(4eg) = 1. For the cluster, we have = 128, i?vr = 20, and fci = 204.8, where Md, 
i?vr, and Eel are the mass, virial radius, and binding energy of the cluster, respectively. The mean velocity dispersion 
of stars in a galaxy is Vg = l/-\/2, and that of galaxies in the cluster is Vc\ = 4/-\/5. The cluster crossing time is 

Ter = ^ = 10^/5~22, (21) 

Vcl 

and the half- mass relaxation time (Spitzer 1987, equation [2-63]) is 

where the cluster half-mass radius i?h = 15.4 and the Coulomb logarithm InA = InA^ = In 128. If we assume that 
= IO^^Mq and r^^ = 30 kpc, the cluster model has = 1.28 x i?vr = 0.6 Mpc, T„ = 2.5 Gyr, and 

tZ = 2.1 Gyr. 

Only models FPAl and FPA2 have a common halo initially (see subsection 5.1); the ratio of the common halo mass 
to the total mass, A/h/Afd, is 1/2 and 3/4, or if the mass is converted into the number of galaxies, iVg = 64 and 32, 
for FPAl and FPA2, respectively. Galaxy and common halo components are distributed in space so that the ratio of 
their densities is independent of the radius. The other conditions are the same as those of the other models. 



4. Simulations 



4-.1. Fokker-Planck Simulations 

Details of the numerical integration scheme of the FP equation have been described by Takahashi (1995, 1997). We 
used 301 energy, 51 angular momentum, and 121 radial grid points. We used K = 20 mass components for galaxies. 
These components are uniformly spaced in m between and 1; mi = l//v, m2 = '2/ K, . . . , niK = 1- In addition, one 
component of toq = was assigned to a background (common halo) . We confirmed that these numbers of grid points 
and mass components are large enough to obtain satisfactory convergence of the results. The relative errors in the 
total mass were less than 0.1% and those in the total energy were less than 0.5% in all runs. 
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We calculated a number of FP models, switching on and off mass stripping and energy dissipation and changing 
the parameters of the stripping rate. These models are summarized in table 1. 

4-2. N-body Simulations 

We compared our FP models with the A''-body models of SFM99 and SFMOO for the same initial conditions. SFM99 
performed a number of simulations for the same initial theoretical model using different random realizations of the 
model and varying the number of particles for representing each galaxy; they confirmed that the results of those 
simulations were almost the same. In section 5, the model PP of SFMOO, denoted by NB in this paper, is compared 
with the FP models. The simulation was carried out using GRAPE-4 (Makino et al. 1997) with the Barnes-Hut tree 
algorithm (Barnes, Hut 1986; Makino 1991; Athanassoula et al. 1998). Each initial galaxy model was represented by 
2048 particles, and thus the total number of particles was 262144. For more details of the numerical integration and 
for a galaxy identification scheme, see SFM99 and SFMOO. 

5. Comparison between Fokker Planck and AT-body Models 

5.1. Models with Common Halos but without Mass Stripping 

First we discuss the evolution of model clusters in which mass stripping from galaxies is not allowed, but which 
initially have massive common halos, in order to reveal how such multicomponent systems evolve owing to two-body 
relaxation only. 

Although coUisional evolution of multicomponent stellar systems has been investigated by many authors using FP 
models in the context of globular cluster evolution (e.g., Inagaki, Wiyanto 1984; Cohn 1985; Chcrnoff, Weinberg 1990; 
Takahashi 1997 and references therein), massless particle components were not considered in those studies. On the 
other hand, it is not very difficult to extrapolate the behavior of massless particles from the behavior of particles 
whose masses are much smaller than the mean particle mass. In that sense, no essentially new results were found 
in our simulations, as discussed below. However, it is useful to see first how "pure" coUisional systems evolve before 
considering more complicated systems with mass stripping. Yepes et al. (1991) and Yepes and Domfnguez-Tenreiro 
(1992) also studied pure coUisional evolution of galaxy clusters including massless particles, but using the moment 
equations of the FP equation. 

Figure la shows the evolution of density profiles for model FPAl; the solid lines and dotted lines represent the 
galaxy and common halo components, respectively. Figure lb shows the corresponding profiles of the logarithmic 
gradient defined by 



Figures Ic and Id are for model FPA2. 

We have assumed that the number of galaxies, iVg, is 64 and 32 for FPAl and FPA2, respectively. However, it 
should be noted that A'g is only concerned with time scaling; if time is measured in units of the relaxation time, the 
evolution of these models docs not depend on the number of galaxies, because only the two-body relaxation process 
is taken into account in them. We set the Coulomb logarithm equal to ln(A/ci/'Tig) = In 128 for time scaling for both 
models. 

As is well known, energy redistribution among particles due to two-body relaxation eventually leads the cluster core 
to gravothermal core collapse (Lynden-Bell, Wood 1968; Lynden-Bell, Eggleton 1980); after the onset of core collapse 
the central density continues to increase and the core radius continues to decrease, unless other physical mechanisms 
intervene. Late core collapse proceeds sclf-similarly with the development of the power-law density cusp, p oc r~°' . 
Cohn (1980) found a = 2.23 for single-component clusters. In multicomponent clusters, the most massive component 
eventually dominates the central part, i.e., strong mass segregation occurs, owing to a tendency toward equipartition 
of energy. Thereafter, the most massive component collapses as in single-component clusters and lighter components 
have shallower density cusps (e.g., Inagaki, Wiyanto 1984; Chernoff, Weinberg 1990). 

The development of the density cusps of the galaxy and common-halo components through gravothermal core 
collapse is clearly shown in figure 1. Power-law cusps appear for r ^1 = i?vr/20. The slopes of the cusps hardly 
depend on the difference in the initial conditions of FPAl and FPA2; the slope of the galaxy cusp approaches a = 2.23 
and the slope of the common-halo cusp is much shallower, a ~ 0.5, as expected. 

Cohn (1985; see also Chcrnoff, Weinberg 1990) found that in deep collapse phases of multimass clusters, the power- 
law indices Ui of the density profiles of different components are approximately related by 



where mu is the particle mass of the most massive component, which dominates the central potential well. This 
equation is obtained on the assumption that the distribution function of component i is given by fi oc {—E)P' with 




(24) 
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Pi/Pu = rai/m^. This relation for pi was originally found by Bahcall and Wolf (1977) for a steady state solution of the 
stellar distribution around a central massive black hole. In deriving equation (24) the results of numerical simulation 
(Cohn 1980), = 2.23 and = 8.2, are also used. 

Equation (24) gives 0.345 as m, 0. Figures lb and Id show that the actual slope of the common halo cusp 
is slightly larger than this predicted value. Such deviation from equation (24) at small masses was also reported by 
Chernoff and Weinberg (1990). 

For models FPAl and FPA2, we followed the cluster evolution up to a rather deep collapse phase in order to see 
self-similar development of the density cusps. However, we should keep in mind that it is very unlikely that actual 
galaxy clusters can reach such deep collapse, because galaxies have appreciable finite sizes compared to the cluster 
sizes. For our model cluster, the sizes of the galaxies are ^ 1, and thus the core radius of the cluster cannot be less 
than ^ 1. 

5.2. Models with Mass Stripping but without Energy Dissipation 

In this section wc show the results of FP models in which mass stripping is included but energy dissipation is 
neglected, and compare these results with iV-body results. 

Figure 2a shows the fraction of the common halo mass to the total mass as a function of time for models FPBl 
and NB. Figure 2b shows the evolution of the central density of the galaxies and that of the common halo, Pg{0) and 
Ph(0). For the A^-body model the galaxy density is not plotted, because it changes very noisily with time owing to the 
small number of galaxies (see figure 6). The center of the cluster of the A^-body model was chosen to be the density 
center (Casertano, Hut 1985) of the common halo particles. 

Concerning the growth of the common halo mass, Mh, the agreement between FPBl and NB is pretty good, though 
the growth speed slows down somewhat earlier in FPBl than in NB. Actually, the value of the numerical constant, 
C = 25, in equation (11) for FPBl was chosen to achieve such a good agreement. The initial growth rate of Mh in 
the FP model is determined by only the stripping rate for given initial conditions. This means that we can always 
adjust C to obtain a desirable initial growth rate. Once its value is chosen, however, we have no freedom to adjust 
later growth of the common halo, which depends on evolving cluster properties. Therefore, the fine tuning of C is not 
responsible for the overall good agreement observed in figure 2a. 

In principle, the value of C can be determined from the data of A^-body experiments of galaxy-galaxy encounters such 
as given by FM. However, its exact value depends on the detailed structure of the galaxies, which changes with time 
under the cluster environment. From equations (8) and (9), we find, for T]y = 2, C" = C(87r/3)-i/Vln(i?vr/rvr) = 0.115C 
if rjr = 2, and C = C{8Tr/3)-^/^/{l/r^, - l/i?vr) = 0.364C if r]r = 3. Hence C = 25 corresponds to C = 2.9 and 9.1 for 
i]r = 2 and 3, respectively. Equation (6) with these values of C is consistent, in order of magnitude, with the data of 
FM (sec their figures 7, 8, 11, and 12). 

Models FPBl and NB arc in rather good agreement in the evolution of A/h, but show a clear difference in pii(O) 
evolution: the central density increases much more slowly in FPBl than in NB. Wc will see below that the inclusion 
of energy dissipation, which makes the treatment of tidal stripping in FP models fully self-consistent, removes this 
difference. 

The initial quick increase in /3h(0) and the initial decrease in Pg(0) are simply due to the fact that mass stripping 
from the galaxies suddenly starts at t = 0. After a few cluster crossing times, the central galaxy density begins to 
increase. Mass stripping works only to decrease the galaxy density by removing part of galaxy mass to the common 
halo. Therefore, the increase in the galaxy density should be due to two-body relaxation. Since the processes of core 
collapse and mass stripping occur simultaneously, the central density of the galaxies and that of the common halo 
stay comparable for a long time, while the central density is soon dominated by galaxies in models FPAl and FPA2 
where there is no mass stripping. Figure 2b shows that the galaxy density finally exceeds the common halo density 
at the center even with mass stripping. We discuss in more detail the process of core collapse with mass stripping in 
section 6. 

Here, we mention the choice of the value of the Coulomb logarithm. As noted above, the central density increase 
in model FPBl is significantly slower than that in model NB. Wc found that rough agreement in the central density 
evolution was obtained when we artificially increased the two-body relaxation rate by about a factor of two without 
changing the mass stripping rate. Increasing the relaxation rate by a factor of two is equivalent to doubling the 
Coulomb logarithm. In A, i.e., replacing In 128 with 2 In 128 = In 128^ in the present case. There is some uncertainty 
for the choice of the exact value of A, but it is very unlikely that such a large value as A = 128^ is justified. 

5.S. Models with Mass Stripping and Energy Dissipation 

We now take account of the effects of both mass stripping and energy dissipation in FP models (sequence FPC 
in table 1). The evolution of model FPCl (our standard model) is shown in figure 3. The growth of Mh in FPCl 
is similar to that in FPBl, although the former shows a slightly larger growth rate at late epochs. The difference 
between FPBl and FPCl is clear in the evolution of the central density, which increases more rapidly in FPCl than 
in FPBl. We now see good agreement between models FPCl and NB in the central density evolution as well as in the 
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evolution of the common halo mass. It should be recalled that energy dissipation is included consistently with mass 
stripping, without introducing any additional free parameters. 

In the following, we investigate how sensitively the behavior of FP models depends on the parameters of the mass 
stripping rate (C, r]^^ rjm, and C). The probable ranges of these parameters are limited by the results of A^-body 
experiments and theoretical considerations (see section 2.2). 

In figure 4 we compare models FPCl, FPC2, and FPC3, which differ only in the value of C (C = 25, 15, 35). Larger 
C gives a higher stripping rate, and consequently faster growth in the mass of the common halo. At early epochs the 
central density of the common halo also increases more quickly for larger C, but at late epochs it increases more slowly 
because faster mass stripping prevents the density increase more severely. We conclude that the value of C = 25 gives 
the best agreement between the FP and A^-body results. 

Figure 5 compares models FPCl, FPC4; FPC5, and FPC6, which differ in the value of ry^, rjm, or ( (see table 1). For 
model FPC4 it is necessary to adjust the value of C so that the initial growth rate of the common halo becomes similar 
to those of the other models. These models are not very different from one another. Therefore, we conclude that the 
results of FP simulations do not very sensitively depend on the mass-stripping formula as long as the parameters are 
in reasonable ranges. 

One noticeable difference is seen in the evolution of the central density of model FPC5 (77™ = 1); at late epochs 
the central density increases more slowly compared with the other models. We find that this is because massive 
galaxies are much more depleted in model FPC5. In all models, the mean galaxy mass {mf)p in the central regions 
first decreases quickly and then stays roughly constant as the galaxy density increases. At late epochs, {mf)p 0.2 
at the center for model FPC5 and ~ 0.4 for the others. This implies that the collapse rate of FPC5 is smaller by 
about a factor of two. Equation (18) tells us why such a difference occurs; tidal stripping time Tts is proportional to 
1/(to/)p for r]ra = 2, but Tts does not depend on {mf)p for rjm = 1- Therefore, while Tts increases as {'mf)p decreases 
in the rjm =2 models, it stays almost constant (actually somewhat decreases) in model FPC5. Consequently, {'mf)p 
decreases to a much lower value in FPC5. A comparison with the A^-body model in figure 5b supports rjm — 2 rather 
than 7y„j = 1. 

Figure 6 shows the mass-density profiles of models NB (left) and FPCl (right) at selected epochs. The density 
profiles for the galaxies (middle) and the common halo (bottom) are plotted separately as well as the total density 
(top). For model NB, the density profiles of the galaxies are rather noisy because of the small number of galaxies. We 
should note that a comparison of the profiles of the TV-body and FP models on a length scale of ^ 1 is, at least for 
the galaxies, meaningless, since the sizes of the galaxies, Tg, is ~ 1. Thus, for r ^ 1 in the iV-body model, we see the 
internal structure of a galaxy if any galaxy is there, and a void if not. At late epochs, the central region of the cluster 
is usually occupied by a few massive galaxies that have sunk there owing to dynamical friction. For r ^ 1, models NB 
and FPCl show very similar evolution of the density profiles, apart from fluctuations in the A'^-body model. 

The similarity between the two models is more clearly shown in figure 7, which directly compares their density 
profiles at < = 8.9 Td-. In particular, almost perfect agreement is seen in the density profile of the common halo. 

The agreement between the FP and A^-body results shown above is enough to convince us that our FP models 
include all essential physical processes concerned and can describe the cluster evolution rather accurately. 

In figure 8 we show the logarithmic density gradient, a [sec equation (23)], as a function of the radius for model 
FPCl at i = 8.9 Ten together with the density profiles of FPCl and NB again. For 1 ^ r < 5, the value of a for the 
total density varies gently with the radius from ^0.5 to ~ 1.5 with a median value of ~ 1. This implies that the total 
density profile in this region may be reasonably approximated by a power law, r~^, as SFM99 claimed. 

The common halo component dominates the density for 1 ^ r ^ 20. For r ^ 1, the density of the galaxy component 
overwhelms that of the halo component in the A^-body model, while they are comparable in model FPCl. As noted 
above, this difference arises because in the A'^-body model a massive galaxy stays around the cluster center at late 
epochs and its internal structure is observed on this length scale. In fact, the density of the particles in each galaxy is 
much larger than the density of the common halo particles for ?' ^ 1 in NB. This extra density increase in the A^-body 
model strengthens our impression that the total density profile is given by a power law. 

It is not clear from figure 8 whether the central density distribution is really described by any single power law, 
since the central density has not yet increased sufficiently from the initial value. In principle, FP calculations can be 
(formally) continued further, if we neglect the fact that the galaxies have finite sizes. However, we could actually not 
continue FP calculations to very advanced stages of core collapse for model FPCl and other similar models, unlike the 
case for the models without stripping, owing to numerical instability. This difficulty of FP calculations is linked to the 
divergence of the distribution function at the center caused by continuous decrease in the central velocity dispersion; 
see section 6. 

In order to see the mass dependence of the galaxy distributions, we divided the galaxies into four mass groups, 
m e (0,0.25], (0.25,0.50], (0.50,0.75], and (0.75,1], and plotted the density profile for each group for model FPCl at 
t = 8.9 Tcr h^ figure 9. It is apparent that heavy galaxies are strongly depleted in the central regions, owing to tidal 
stripping, and have a very flat density proflle. If there were no stripping, the heavy galaxies would dominate the 
central regions as the result of two-body relaxation. 
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The development of "inverse mass segregation" is clearly shown in figure 10, which plots the mass- weighted mean 
galaxy mass, {m)p, as a function of the radius for model FPCl. The mean mass decreases toward the center. This 
indicates that mass stripping is more efficient than dynamical friction for heavy galaxies. Only at very late epochs 
(t ^ 8.9 Tcr), when the central density of the galaxies has increased sufficiently high, is a slight increase in the mean 
mass seen at the center. Similar inverse mass segregation was already found by Funato et al. (1993) and SFM99 in 
A''-body simulations and by Merritt (1984) in FP simulations. 

The density profiles of models FPC1-C6 are similar, if they are compared at epochs with the same central density. 
In particular, the halo density profiles are almost indistinguishable. However, we show in the next section that the 
density cusp gradient does depend on the stripping rate. It is obvious that the gradient, a, should be close to 2.2 
when mass stripping is almost negligible compared to two-body relaxation. On the other hand, if stripping is infinitely 
fast, all galaxies will dissolve instantaneously and no further evolution of the cluster can occur based on our present 
assumptions. The stripping rates for models FPC1-C6 are not very different (see figures 4 and 5), and hence the 
difference in their density profiles is insignificant. 

6. Cool Core Collapse and Shallow Density Cusps 

In subsection 5.3 we showed that a shallow density cusp, approximated by p oc r~^, develops in the central regions 
of a galaxy cluster. Now let us assume that a collapsing core leaves behind a cusp that exactly follows a power-law, 
/OCX r~". This implies that the central density grows as pdt) oc r~"(t), as the core radius rdt) decreases with time. 
Similarly, a dimensional analysis implies that the central velocity dispersion changes as v'^{t) oc Gpcr^/rc oc r^~"(t) 
(but this is valid only for a > 1 as we discuss later). This relation predicts that the central velocity dispersion increases 
(i.e. the core gets hotter) as the core shrinks if a > 2, and that it decreases (i.e. the core gets cooler) if a < 2. For 
core collapse in star clusters, detailed numerical calculations (Cohn 1980; Lynden-Bell, Eggleton 1980) show that 
a = 2.2 > 2, and that actually both the central density and velocity dispersion increase as the core shrinks — "hot 
collapse" occurs. For core collapse in galaxy clusters, our calculations indicate a < 2. Therefore, we expect "cool 
collapse" in this case. 

Lynden-Bell and Eggleton (1980) proved that a should be in the range 2 < a < 2.5 for the existence of any self- 
similar collapse solution in equal-mass clusters. The condition a > 2 comes from the requirement that temperature 
must decrease outward so that the core can shrink by losing heat to the halo. In multimass clusters, however, heat 
is also transferred among different mass components. Therefore, it can happen that the core of heavy particles loses 
energy to light particles and shrinks even if a < 2. 

Let us see in more detail the case of a shallow density cusp, 

p{r)(xr~" (0<a<2). (25) 

This produces a potential, 0(r), that is finite at the center and satisfies 

</>(r)-0(O)ocr2"". (26) 

If we assume that the velocity distribution is isotropic (this is not a bad assumption near the center), we find that the 
density distribution (25) corresponds to the distribution function 

/(i;)oc [i;-0(O)]^(^~")/2(2-a)^ (27) 

Thus, / diverges as i? ^ '/'(O). This divergence causes difficulty in numerical integration of the FP equation, and 
therefore we had to stop our calculations at relatively early (low central density) collapse phases for the models with 
mass stripping. The depth of the central potential well becomes shallower as a decreases (Tremaine et al. 1994). 

As is argued by Tremaine et al. (1994), the dependence of the velocity dispersion, w^^, on the radius for the density 
cusp represented by equation (25) changes at a = 1. We find 

i;f„(r) ocr^-", for a > 1, (28) 

which is the equation expected from a dimensional analysis. For a < I, the behavior of the velocity distribution is 
more subtle, since it is dominated by high-energy parts of the distribution function, which must deviate from equation 
(27) if the total mass is to be finite. If we consider any finite-mass model, we find 

?;^(r)ocr" for a < 1 (29) 

(see also Dehnen 1993). For any a < 2, ^ as r — > 0, and hence cool core collapse is expected. 

In figure 11 we plot the evolution of the central densities and the central velocity dispersions of the galaxy and 
common halo components for models FPAl and FPCl. Here the velocity dispersion of the galaxies is the mass- 
weighted average over all the galaxies. The central density increases with time in both models, but the behavior of 
the velocity dispersion is very different between FPAl and FPCl. As we expected, the central velocity dispersion 
of the galaxies does decrease as core collapse proceeds in model FPCl, while it increases in model FPAl. In model 



Evolution of Galaxy Clusters 



11 



FPCl, even the central velocity dispersion of the common halo finally decreases when core collapse accelerates. It 
is considered that this is because cool halo particles are being supplied from cool central galaxies. (Remember we 
assumed that the particles stripped from a galaxy initially have the same energy and angular momentum per unit 
mass as those of its parent galaxy.) 

We now consider why the velocity dispersion of the galaxies continue to decrease when mass stripping is included. 
Spitzer (1969) showed that equipartition among different mass components cannot be achieved if self-gravity of heavy 
components is strong. It is instructive to recall here his simple model. Consider a spherical system composed of 
particles of two different masses, mi and m2 {mi < 7712). We assume that the total mass of the light component is 
much larger than that of the heavy component, Mi 3> M2, and that the heavy component is concentrated at the 
center. Then, the virial theorem for component 2 can be written as 

2 GM2 , ^ GMiRl 

-2=a^ + b^^, (30) 

where V2 is the mean velocity dispersion for component 2, Ri and R2 arc the median radii of components 1 and 2, 
respectively, and a and b are numerical constants. For fixed Mi, M2, and Ri, the right side has a minimum value of 

2 ^1 2.M/3 G(A/iA/|)V3 

'^2.min = 2(a h) ' (31) 

Hi 

at R2 — R2,cr\t = 

{aM2lhMif/^Ri. The minimum of z;| exists because the heavy particles form a self-gravitating 
system by themselves when R2 i?2,crit, but behave as a non-self-gravitating system when R2 S> i?2,crit- Therefore, 
after the concentration of the heavy component becomes sufficiently high, its velocity dispersion increases as the heavy 
component loses energy and contracts. However, if we now allow M2 to change with Mi + M2 fixed, we find that 
V2 niin and i?2,crit ^ as M2 0. This implies that the velocity dispersion of the heavy component can continue 
to decrease to zero as the core shrinks, if mass stripping proceeds at an appropriate speed. Mass stripping prevents 
the heavy component from forming an independent self-gravitating system. 

We may expect that core collapse changes its nature from cool collapse to hot collapse as tidal stripping becomes 
less and less important compared with two-body relaxation. To confirm this, we ran models FPC7, C8, and C9, which 
are the same as model FPCl but have C = 2.5, 0.25, and 0.025, respectively. In principle, we should change the initial 
models of galaxies and clusters rather than C, which is considered to be roughly constant, in order to see the evolution 
of clusters with different stripping importance. From equation (20) we find that stripping becomes less important as 
Wgo/Vci increases with the other parameters fixed. Thus, very low stripping rates imply the condition Vgo 3> Vd- (cf. 
In star clusters the escape velocity from a star is much larger than the velocity dispersion of stars.) However, we 
should not directly apply the equations for tidal stripping presented in subsections 2.2 and 2.3 to such cases, since 
those equations are considered to be valid only for Vd ^ Vgo- We therefore choose to change C just to conveniently 
control the stripping efficiency in our theoretical experiments. 

Figure 12 shows the growth of the common halo mass for models FPCl and FPC7~C9. It is interesting that the differ- 
ence in the growth rates at the late times is within an order of magnitude among these models {[Tcr/Md] [dMi^/dt] ^^0.03, 
0.08, 0.05, and 0.01 for FPCl, C7, C8, and C9, respectively), while C varies over three orders of magnitude. This may 
be explained by a kind of self-regulation mechanism: if C is large, the galaxies quickly become more compact through 
initial violent stripping, and then the stripping speed slows down; if C is small, the central density can increase to 
very high values, and this finally increases the stripping rate at the central regions. 

Figure 13 shows the evolution of the central density and velocity dispersion of the galaxy component for models 
FPCl and FPC7-C9. The transition from cool collapse to hot collapse occurs at C ~ 0.25. In model FPC8 the central 
velocity dispersion increases until the central density increases by about a factor of 10"^, but decreases after that. 
Equation (20) tells us that in the present models the stripping time eventually becomes shorter than the relaxation 
time as Vd increases toward infinity. Thus, we expect that for any small value of C, tidal stripping will become 
effective and the velocity dispersion will begin to decrease when the velocity dispersion has increased to a critical 
value, which is inversely proportional to C. The above results show that, unless the stripping rate is decreased to a 
very low value compared to the standard value, cool collapse occurs. Therefore, under the conditions of real galaxy 
clusters (Vci ^ "^^go); cool collapse is expected to occur. 

Figures 14 shows the profiles of the density and the logarithmic density gradient a for model FPC7 (C = 2.5), and 
figure 15 shows those for model FPC9 (C = 0.025). In model FPC7 a shallow density cusp of a < 2 develops, but 
this cusp is steeper than the cusp in model FPCl (see figure 8). In model FPC9 a deep cusp of a > 2 develops, and 
a approaches 2.2 as the core collapse proceeds. Unlike in no-stripping models FPA, the central density of the halo 
component also significantly rises in model FPC9 due to continuous production of new halo particles. 

Figures 16a and b show the evolution of velocity dispersion profiles for models FPC7 and FPC9, respectively. In 
these figures qualitative difference between models FPC7 and FPC9 is clearly seen. In model FPC7 an inverse velocity 
dispersion gradient develops in the central regions. 
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7. Summary and Discussion 

In this study, wc investigated the dynamical evolution of galaxy clusters after virialization using FP models and 
compared these models with the A^-body models of SFM99 and SFMOO. Our FP models included the effects of 
gravitational two-body encounters between galaxies and between galaxies and common halo particles which are much 
lighter than the galaxies. Furthermore, tidal stripping of mass from the galaxies to the common halo and accompanying 
dissipation of the orbital kinetic energies of the galaxies were included using revised cross-section formulae. 

We found that the results of the FP models agree very well with those of the A^-body models regarding the growth 
of the common halo mass and the evolution of the cluster density profiles. This indicates that all of the important 
physical processes that occur in the A^-body models are properly included in the FP models. Although there is some 
ambiguity in the cross sections of mass stripping and energy dissipation, we found that the results of the FP simulations 
do not sensitively depend on the details of the cross-section formulae. 

In the evolution of galaxy clusters, tidal mass stripping from galaxies is very important. Under the initial conditions 
that SFM99 and SFMOO employed, where the total cluster mass is initially attached to galaxies, more than half of the 
total mass turns into the common halo during a first few crossing times. Heavy galaxies sink toward the cluster center 
owing to dynamical friction. However, a high density of galaxies results in strong stripping, and thus it is usually 
seen in our models that the mean mass of galaxies decreases toward the center, which is inverse mass segregation. 
Dissipation of the orbital kinetic energies of galaxies due to tidal stripping also plays an important role in accelerating 
the sinking of the galaxies toward the cluster center, if the velocity dispersion of the stars in the galaxies is not very 
small compared to that of the galaxies in the cluster. 

As the central density increases, a cusp profile, approximated by a power law, pcx r~", develops at the core region. 
In the model clusters of SFM99 and SFMOO, a for the total density is This density distribution consists of the 
galaxy component with larger a (still < 2) and the common halo component with smaller a. The development of the 
density cusp is a consequence of the collisional evolution of galaxy clusters. 

The slope of the density cusp, a, depends on the ratio of the stripping rate to the relaxation rate. When the effect of 
stripping is negligible compared to relaxation, "hot core collapse" , which is usual gravothcrmal collapse, occurs: i.e., as 
the core shrinks both the central density and velocity dispersion increase, and a deep cusp with a ~ 2.2 > 2 develops. 
When stripping is important, "cool core collapse" occurs: while the central density increases, the central velocity 
dispersion decreases, and consequently a shallow cusp {a < 2) develops. Faster tidal stripping results in shallower 
cusps. Under the conditions of real galaxy clusters, tidal stripping is rather effective and cool collapse is expected to 
occur. 

One might think that this conclusion is in contradiction with the results of SFMOO, who showed using A^-body 
simulations that the profile with a ^ 1.2 (1-1.5) develops regardless of the initial models of galaxies and clusters. 
Actually, their results arc consistent with the results of this paper. SFMOO varied the initial models of galaxies and 
clusters adopting Plummer, King (King 1966), and Hernquist models (Hernquist 1990), but did not vary the number 
of galaxies and the ratio of the galaxy virial radius to the cluster virial radius. As a result, the growth rates of the 
common halo mass, or the mean stripping rates are similar among all of their runs. Our FP simulations show that 
while a depends on the stripping rate, it does so only weakly. Therefore, it is no wonder that no significant difference 
in a has been detected among the A^-body models of SFMOO. 

Hot core collapse terminates with the formation of hard binaries. Because of subsequent energy release from these 
binaries, the core starts to expand (see, e.g., Spitzer 1987). What is the end state of cool collapse? Through cool 
collapse, low-energy galaxies accumulate around the cluster center. They finally merge and form a massive central 
galaxy. The central galaxy will further grow by eating galaxies falling to the center in succession, and this process will 
continue until all galaxies, except for the central one, disappear. Thus, we will see at last one giant (cD ?) galaxy with 
an enormous halo. In the present study we did not include the merging process in FP models, and hence we could 
not follow very late stages of the cluster evolution with FP models. Future A^-body calculations as well as improved 
FP models will reveal the consequences of cool collapse in more detail. 

This work was supported in part by the Research for the Future Program of Japan Society for the Promotion of 
Science (JSPS-RFTP97P01102). 
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Fig. 1. (a) Evolution of the density profiles of the galaxy component (solid) and common halo component (dotted) for model FPAl. 
Plotted are the profiles at t/Tcv = 0, 6.7, 7.99, 8.0209, 8.0214. The central densities of both components increase with time, except 
that the central density of the halo component decreases initially, (b) The logarithmic density gradient cx = —dlnp/dlnr for the 
density profiles shown in (a). The two horizontal dashed lines are a = 2.23 and a = 0.345 (see text), (c) Same as (a), but for model 
FPA2. Plotted are the profiles at t/Tcr = 0, 4.5, 5.42, 5.4454, 5.4457. The central densities of both components increase with time, 
(d) Same as (b), but for model FPA2. 
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Fig. 2. (a) Ratio of the common halo mass, A/jj, to the total 
cluster mass, Md, as a function of time. The circles repre- 
sent model NB and the solid line represents model FPBl. (b) 
Evolution of the central density of the common halo compo- 
nent (solid) and that of the galaxy component (dashed) for 
model FPBl. The circles represent the central density of the 
common halo for model NB. 



0.8 

i^"" 0.4 
0.2 





1 1 1 

- (a) 


1 1 1 1 1 


' 1 ' 


- 






- /o 






_ /o 






/o 






h 

^111 


1 1 1 1 1 


, 1 , 



5 10 

-1 1 1 1 1 1 1 1 1 1 r- 

(b) 




Fig. 3. Same as figure 2, but model FPCl is compared with 
model NB. 
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Fig. 4. Same as figure 2, but the results for models FPCl, 
FPC2, and FPC3 arc shown by the solid, dashed, and dot- 
ted lines, respectively. In (b) only the central density of the 
common halo is shown. 



Fig. 5. Same as figure 4, but the results for models FPCl, 
FPC4, FPC5, and FPC6 are shown by the solid, dashed, dot- 
ted, and dash-dotted lines, respectively. 
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Fig. 6. Density profiles of models NB (left) and FPCl (right) at selected epochs. The top panels show the total density, and the 
middle and bottom panels show the densities of the galaxy component and the common halo component, respectively. The density 
profile of the initial cluster as well as the profile of an initial galaxy (here its center is placed at the cluster center) is plotted with 
thin dotted lines. 
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Fig. 7. Comparison of the density profiles of models NB 
(solid) and FPCl (dotted) at t = 8.9Tcr. The top panel shows 
the total density, and the middle and bottom panels show 
the densities of the galaxy component and the common halo 
component, respectively. Note that the two models should be 
compared with caution for r ^ 1, because in model NB the 
galaxies have sizes of ~ 1 and their internal structure is seen 
on the scale of r "5 1 . 
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Fig. 8. (a) Density profile p(r) for model FPCl at t = 8.9Tcr. 
The solid, dashed, and dotted lines represent the total den- 
sity, the galaxy component, and the common halo component, 
respectively. The density profiles for model NB at t = 8.9Tci- 
are also shown by thin lines, (b) Logarithmic density gradi- 
ent, a{r) = —dlnp/dlnr, corresponding to the density profiles 
for model FPCl shown in panel (a). 
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Fig. 9. Density profiles for the galaxies in four mass ranges, 
(0,0.25], (0.25,0.50], (0.50,0.75], and (0.75,1], separately 
shown for model FPCl at t = 8.9Tcr. The dotted hne is for 
the common halo. 
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Fig. 11. (a) Evolution of the central densities of the galaxy 
component (dashed) and the halo component (dotted) for 
models FPAl and FPCl. (b) Evolution of the central ve- 
locity dispersions of the galaxy component (dashed) and the 
halo component (dotted) for models FPAl and FPCl. 



100 



Fig. 10. Mean galaxy mass {m)p as a function of radius for 
model FPCl at selected epochs. 
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Fig. 12. Growth of the common halo mass for models FPCl, 
FPC7, FPC8, and FPC9. 



10' 



100 



o 



0.1 
0.01 



10 



-3 



= 1 1 1 1 1 1 1 1 

= (a) 


1 1 1 1 


I'll 

: 1 1 


1 FPn 






L _ _ FPC7 




' ' 1 


! FPC8 






;. FPC9 


/ 






















Illlllll 


1 1 1 1 


III" 







5 10 
t/T 



15 



cr 



10 



o 

c\l 



- lllllll 

-(b) 


1 1 1 1 1 1 1 - 


Illlllll 


\ : 

V 1 

, 1 , , , , 1 " 







5 10 



15 



0.1 F 



0.01 



10' 



: 1 1 1 1 1 1 ll| 1 1 1 1 1 1 1 1| 

'- ^^'^ 


1 1 1 1 1 1 ii| 1 lllllll 

FPC7 ^ 








1 1 





-3 
-4 
-5 

0.01 0.1 



1 

r 



10 100 



2 - 



0.01 




10 100 



Fig. 14. (a) Evolution of the density profile for model FPC7. 
The profiles at t/Tcr = 7.1, 9.8, and 10.3 are plotted. The 
solid, dashed, and dotted lines represent the total density, the 
galaxy component, and the common halo component, respec- 
tively, (b) The profiles of the logarithmic density gradient, a, 
corresponding to the density profiles shown in panel (a). 



Fig. 13. (a) Evolution of the central density of the galaxy 
component for models FPCl, FPC7, FPC8, and FPC9. (b) 
Evolution of the central velocity dispersion of the galaxy com- 
ponent for models FPCl, FPC7, FPC8, and FPC9. 
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Fig. 15. Same as figure 14, but for model FPC9. The profiles 
at tjTcT = 7.8, 12.7, 14.3, and 14.7 are plotted. 



Fig. 16. (a) Evolution of the velocity dispersion profiles for 
model FPC7. The profiles at the same epochs as in figure 
14 are plotted. The dashed and dotted lines represent the 
galaxy and common halo components, respectively. For the 
galaxy the mass-weighted mean velocity dispersion is plotted. 
The central velocity dispersions of both components decrease 
with time, (b) Same as (a), but for model FPC9. The profiles 
at the same epochs as in figure 15 are plotted. The central 
velocity dispersions of both components increase with time. 



